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.SECTION I- INTRODUCTION 

This final report on Task F of Contract NAS6-2307 must be regarded as 
a report of progress to date in the problem of analysis and Interpretation 
of the waveform-related data from the S-193 radar altimeter experiment. 

The waveform analyses are not complete and a continuing series of different 
”cut and try” calculations must be done; the procedures which are presented 
In this report provide the basic building blocks for this continuing search. 

In very brief summary of the results, two major problem areas exist. 

1.) determining the sample-and-hold <S&H) dc offsets, and also the parameters 
for a theoretical radar mean return waveform from a set of averaged S&H 
data as input, and 2.) performing a time-realignment of individual sets of 
SSH results prior to averaging in order to correct for the altimeter tracker 
time-jitter. For the first problem area, a general-purpose chi-squared 
minimizing computer program, X2MIN, has been used- to best— fit the non- time- 
realigned S&H averages, and the results are encouraging - a set of S&H offsets 
Is found and some consistency is observed between these offsets seen in 
various ”good'* data passes. Time-realignment and averaging, the second prob- 
lem area, is not in such good shape; we have thus far been unable to produce 
any improvement In our knowledge of the waveform by use of time-realignment 
averaging, and therefore we have no confidence that any need exists 
for applying this cumbersome procedure to waveform data. These two areas, 
line-fitting and time- realignment, will be more thoroughly defined and dis- 
cussed in the subsequent sections of this report. 

For the waveform experiments the altimeter’s tracker must be operating 
("locked”) when S&H data are acquired. This eliminates from any further con- 
sideration here the dual-pulse mode. Mode 3, and the Nadir Align Mode. Table 1 
summarizes the (nominal) pulsewldths, IF bandwidths, S&H positions and spacing, 
and the number of frames of data for the remaining altimeter modes and submodes 
for which the altitude tracker is locked. 

Virtually all our effort for this report has been spent on Mode 5, 
subraode 2, the "brute force” short pulse (a nominal 20 ns pulse) experiment 
which was intended as a direct comparison submode for Mode 5’s pulse compres- 
sion submode, SMl. Since M5,SM1 was apparently operating Incorrectly until 
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late In mission SL-3, only M5,SM2 remain as a source of short-pulse ocean- 
scattering information. The longer-pulse submodes (100 ns nominal) are im- 
portant for determination of the ocean’s radar backscattering cross-section, 
Cf*^* whose effects are more pronounced in the trailing edge of the mean return 
waveform; apart from noting that the waveform fitting procedures described 
later in this report may also be used for refining the antenna pointing angle 
estimates which are necessary for calculations, we will not deal further 
here with a°. If any changes in rms ocean surface roughness are detectable 
in the S-193 altimeter data, these changes should manifest themselves in 
the leading edge of the short-pulse mean return waveform. Moreover, the 
tracker jitter should leave no measurable effect on the leading edge of the 
100 ns mean return waveform while estimates made prior to Skylab launch indi- 
cated that the tracker jitter would have an appreciable effect (increasing 
the apparent risetlme and also increasing the variance) on the leading edge 
of the 10 ns short pulse mean return waveform. (Note the distinction here 
between the nominal 10 ns design goal and the nominal 20 ns actually realized 


pulses; the tracker jitter corrections - the time-realignment procedures - 


Figure 1 sketches the overall waveform processing carried out by a 
series of labelled boxes; this is the diagram appropriate to an automatic 
waveform program which does not yet exist. Because of the number of dif- 
ferent difficulties with the S-193 output data and of the only mixed successes 
in our work with the best of this data, such an overall program not only 
doesn’t now exist but never will. 

For Instance, Box A of Figure 1 refers to editing of data, but because 
of the number of missing modes or submodes, of tracker loss-of— lock, and 
of other data drop-outs, it has not been possible to even begin to define an 
automatic data editing set of criteria. Instead the editing has been a manual 
operation; in fact only those submodes having no apparent loss-of-lock or 
data dropouts over an entire submode have been used in our work to date, 
deferring until later the questions of what to do about cases in which part 
of a submode should be edited out. 

The waveform time— realignment procedure includes the functions of boxes 
E,F, and G of Figure 1 to produce the necessary tracker jitter time-correction. 
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Figure 1. Idealized Overall Waveform Processing Procedure 
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then boxes H and J use these corrections to obtain the final time-realigned 
average waveform. At present, these two steps (that is, D&E&F and then G&J) 
are performed separately instead of within a single program. Box K is an 
eaitirely straightforward application of the line fitting program to the final 
time-realigned average waveform sample but has not been implemented because 
of difficulties to be discussed in the section on time-realignment. The 
programming for box L of Figure 1 has been completed and is sketched later 
In this report; at the time of this report's preparation the histograraming 
had not yet been done on altimeter waveform data. 

The two following sections of this report describe first the line- 
fitting process, then thtitime realignment process. The computer programs 
developed or modified for these tasks are described briefly in the appropri- 
ate sections but the bulk of the computer program documentation is relegated 
to Appendix sections at the end of this report - source decks for all of these 
programs have been separately provided to Wallops Flight Center as these 
procedures have developetiV: through the contract period. In some cases (notably 
the time-realignment main- programs) further program modifications have been 
made after the programs have been installed at the Wallops computer; thus 
the programs now at Wallops vary somewhat from those in the Appendix sections 
of this report. 
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SECTION II. DETERMINATION OF S&H OFFSETS AND LINESHAPE PARAMETERS 

The overall waveform processing diagram, Figure 1 of Section I, shows 
two different uses for the line fitting procedure to be described in this 
section; first, the submode S&H averages are used to find the Individual 
S&H dc offsets, and later the final time-realigned averages are to be used 
to find the final set of descriptive parameters to characterize the waveform. 
The same general-purpose chi-squared-minimizing computer program, X2MIN, is 
used for both of these purposes. The program X2MIN is described in some 
detail in Appendix A, and Appendix B provides, the function subroutine FX 
which is used with X2MIN to fit S-193 S&H data to the theoretical mean re- 
turn waveform described later in this section. 

First we will describe the S&H offset problem in general, and point 
out the lack of a good means to determine the needed S&H offsets from any 
of the calibration modes. Then after describing in general the method to 
be used, we will discuss the particular functional form which represents 
the theoretical mean return waveform. Following the pages on waveform fit- 
ting to the 4-parameter function plus S&u offsets, we present examples of 
vaveforms and offsets determined in this way. Finally a means of examining 
the statistics at one or several points on the mean waveform is described. 

Because the video output in the S— 193 radar altimeter is ac— coupled to 
the S&H gates, the average of a set of readings from a S&H gate which is 
ggiPpT a noise— only region of the mean waveform (i.e., a point prior to 
the leading edge of the waveform) should be zero. Because the S&H gates are 
less than ideal however, the average is non-zero In practice for a S&H gate 
with a noise— only input. Each S&H gate will differ from others, with the 
result that each S&H gate will have its own non-zero average output for 
noise-only input. Furthermore, the offset for each S&H gate may be dif- 

for changes in any of the followings IF bandwidth, receiver tempera- 
ture, and S&H gate width. It is Important to note that the dc-offset for 
each gate must be subtracted from all values obtained by that particular gate 
before any time-realignment and waveform averaging can be done. [As we will 
see In the section on time-realignment, any partlcxdar specified point, or 
”bln", on the mean return waveform receives contributions from not one but 
several different S&H gates as they are carried past the point of Interest 
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by the tracker jitter.] 

In this work, we assume that the offset characterizing each of the 
8 S&H gates is a simple additive effect and that once the offset is deter- 
mined, all the S&H data may be corrected by simply subtracting the appropri- 
ate offset from each S&H data point. This necessary assumption could very 
well be Incorrect; for example, the offsets might depend somewhat upon the 
level of the signal in, or might exhibit some sort of hysteresis-like effects 
or could have other, more complicated behavior. Any case but that of the 
simple additive offset becomes virtually impossible to analyze from the 
limited amount of data available. 


To determine the S&H offsets under the simple additive assumption, we 
should look at the output for a S&H array position in which all eight gates 
are sampling noise. Immediately we are in trouble for submode 1 and 2 of 
Mode 5, as there is no subsubmode satisfying this noise-only requirement. 


The next possibility would be to use the calibration data step (CDS) 
data. Mode 5, submode 3 provide the CDS data for the short pulse mode. Mode 5. 
submode 2, and Figure 2 shows the results given by this CDb ..w.. w.. - 


Pass 9. The solid line in the curve is a Gaussian which is best-fitted to 
the CDS data points. Here again, there la no subsubmode in the CDS data 
in which all eight S&H gates are in a non-signal region. Another diffi- 
culty with all the CDS data in our experience is that it is much less vari- 
able, much more reproducible pass-to-pass than is any of the actual data 


acquisition step (DAS) data In submode 2 of Mode 5. Figure 2 does, however, 
suggest an approach which might be used if we were to use CDS data to deter- 
mine offsets; we could take the offsets as being the differences between the 
actual data points and a best-fitted theoretical curve. Figure 3 shows 
the CDS data of Figure 2 after subtracting the offsets determined in this 
manner, together with a new fitted Gaussian. The fit is somewhat 
better In Figure 3 than in Figure 2 . 


This reasoning can be extended to DAS data as well- All our expecta 
tlons are that the plateau region of the mean return waveform should be a 
smoothly varying function; consequently variations of S&H data about a smoothly 
varying fitted function might be attributable to S&H offset. The functional 
form to be fitted is based on earlier work at Applied Science Associates, and 
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wlll be described in the next paragraphs, but there remains the question of 
where to derive the offsets. Should we take simply subsubmode 3, the set 
of S&H positions starting latest in time on the M5, SM2 waveform? Or use 
some weighted combination of offsets derived from several subsubmodes? The 
answer to which we have evolved is that the eight S&H offsets should be 
treated as free parameters to be varied together with the four lineshape- 
related parameters (A, t^, 3, and C in tlie following paragraphs) for a 12- 
parameter fit to the 32 data points from the 4 subsubmodes of M5, SM2. This 
might seem to overlook the entire subject of time-realignment, since we are 
simply averaging all S&H data over the entire submode but since the primary 
effect of tracker jitter on the observed mean waveform should be to increase 
the leading edge rise time (3 in the 12-parameter fit) , this should not be 
a problem in offset determination. 


The 32 waveform averages y. and their corresponding times t . are as— 
signed an index j=l,32 according to the following rules: j=l,2,...,8 is 

assigned to the result from S&H gate 1,2,..., 8 in sub-submode 0; j«9,10, . . - ,16 
Is assigned to the sub-submode 1 result from S&H gate 1,2,. ..,8; and this 
sequential numbering is continued for sub-submodes 2 and 3. Thus arranged, 
data point j,Cj=l*32), will have come from S&H gate k,(k-l,8), where the re- 
lationship between j and k is k=l + j (mod 8) . The theoretical waveform will 
be a function of four parameters (A,t^,3*C) to be described in the following 
paragraph, and of the eight S&H offsets 0^, (k=l,8). 


The four— parameter waveshape function of time t, f (t) , is derived from 
an earlier report [L. S. Miller and G. S. Brown, "Engineering Studies Re- 
lated to the GEOS-C Radar Altimeter, Final Report for Task D, "NASA Contract 
NAS6-2307, Applied Science Associates, Inc., Apex, N. C, , May 1974, pg. 70], 
and is given by 

£(t) = AP ( g- ° -) exp [^(t - t^)cos|2Cl] sin|25l) , 


for t > t , and by 
o 


t - t 

f(t) ” AP ( g - °) . 


for t < t . In these expressions P(z) is the probability function, which 
— o 

Is related to the error function erf(z) by 
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P(z) - 1/2[1 + erf(-^) ] , 

/2 

and I (z) is the Bessel function of the second kind, order zero. The other 
o 

symbols have the following meaning: 
t * time, in ns 

A « a general amplitude factor, 

t " a reference time origin In nanoseconds 

o 

^ “ antenna angle off-nadir, 

Y “ antenna beamwidth-related constant 

» .00055 (ns)~^^^ for the S-193 radar altimeter (for 

SL-2 and SL-3; the beamwidth apparently changed during SL-4) , 
6 * 2v^c7h 

- 1.656x10”^ (ns) ~ for S-193 

3 “ a pulsewidth-risetime parameter 

The parameter 3 is the product of two parameters CT^ and a of the earlier 
report, where 

a - /Tm PW 

P 

PW « the 3dB pulsewidth of the S-193 altimeter, 
c * the speed of light 
“ .3 m/ns 

and O - rms ocean surface roughness. 

5 

Thus for a flat sea, 3 will be equal to = 7.66 ns (for a 18 ns pulsewidth 
as in submode 2 of Mode 5), and an increase In sea surface rms roughness will 

increase 3 * . 

The jth data point, for the data arranged by submode in the manner al- 
ready described, is where we let y^ be the averaged waveform sample. 

The theoretical function to be fitted will be 

W ‘ ^ ' 

Here the subscript o denotes "observed”, f(t^) is the four-parameter function 
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of A,t *3*C given above, and 0. is the offset of S&H date k with k and j 

related as before. We use a chi-squared minimizing computer routine X2MIN to 

best-fit f to the (y. ,x,). This computer subroutine is described in 
o j j 

Appendix A ; X2MIN requires a function subroutine FX which evaluates the 

function f at each of the input times t.. FX also evaluates the derivatives 
o J 

of f with respect to each of the fitting parameters; these derivatives are 

o 

evaluated at each of the 32 input times t^. 

We therefore need the derivatives of f^ with respect to A,t^,B,C» and 

to 6 . k = 1,8 at each point t.. To handle the offsets, first we use the 
k J 

relationship between k and j to write 

df (t J =1 , if k = 1 + j(mod 8) 

To° J 

k 

"0 otherwise 

Then since f (t^) and f(t.) are related by a simple additive constant (0.), 
o j j 

d_f (t J = df (t ) 

Tp“ 1 dp J . 

where p is any one of the parameters A,t^,B,^. Doing the differentiation, 
and summarizing here, 

- AP(-i^)exp[::|-(tjt^)cos|25li °k • yS 






d,f (t.) « 

K ‘ 


-A 


2tt B 


exp [ ' — 5 1 

2B 


- M ^ 


I 
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A expl^ cos|25| (t.-tj,)] { P(^^) I Ificos|2Cjl^(..-.) 


- sin|25i , , ,, ^ ,l 

/t-=t- ^ 2g2 J/' 


t. > t 

i <> 


At tf \ "A(t -t ) 

— expf 

if - ^ p2 


26^ ’ ■ 


,,t -t . (t.-t )‘ ^2 , 

e [- .1 ^ ^]exp r:f (t -t )cos|2?ll I (...) 
6^ If Y J o ® 


and dfg(t ) = 0 . t < 

il 

• AP(-^j^) '^V%” [“y (tj-t^)cos|25l Jsgn(S) I 2cosl25llj^(...) 


+ S.'t^ sin|2e|l^(...) } . tj > 


In the above, sgn(^) is the sign function, 

8gn(Q = +1 , C ^ 0 

sgn(C) = “I 5 < 0 . 

and I (...) and I, (...) have the same argument (— ✓'t .-t_ sin|2^|) in all 
O 1 T j o 

the above expressions. 

The function subroutine FX is written for the 12 parameters to be stored 

in an array A(I) with the order of parameter assignment as follows: 

A(l) •*- A, amplitude, 

A(2) ^ t , time origin, 

A(3) ^ 3? rise time 

A(4) ^ , antenna off -nadir angle, 

A(5) ^ 0 , offset of S&H //I 

A(6) ^ Oj offset of S&H #2 

• • 

A(12) Oq , offset of S&H #8 

o 
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Appendix B provides a source listing for FX: a flow diagram is also 
provided to simplify program debugging or modification. Source listings 
are also provided for subroutines 5ESI to calculate l^C*) 
for NDTR to calculate P(.). FX includes the option of subtraction of the 
Individual gate offset from the input waveform data after finding those 
offsets. This option is enabled by INFLAG(4) - 1, and assumes the use of 
the REPEAT and the continuation flag on C (refer to Appendix A) , and it is 
assumed that the following data deck organization will be employed. 


Input 

Deck 

1 


1 Comment card 

1 Size card - NP=32,NX=1,NA==4,NC^2, default CHILIM and MAXITR 

1 Data Labels - ’TIME*, ‘WAVEFORM*, * SIGMA* 

32 Data Point cards - (32 input data points) 

4 Variable Par am. cards (A) -initial neighborhood guesses for A(I) ,1=1,4 

2 Const. param(C) - 6,y 


2nd 


inpuc i 

Deck 


[; 

i ^ 

L 


Comment 

Size card - NP=32,NX=1,NA=12,WC=2, default CHILIM and MAXITR 
Data Label Card - REPE in columns 1-4 
12 Variable Param (A) - eontinue Flags on, A(l) , . .. ,A(4) ,zero a.i 

first guess for A(5) , . . .,A(12) 

2 Const. Param. (C) - contunue flags 


3rd 

Input 

Deck 


1 Comment 

1 Size card-NP=32,NX=l,NA=4,NC=10, default CHILIM and MAXITR 
1 Data Label - REPE in columns 1-4 
4 Variable Param. (A) - continue flags 
10 Const. Param. (C) - continue flags on all zeros 


With this organization of the data deck for any one set of input averaged 
waveform points, there will be three successive sets of problem output and 
these will provide the following: 


Prints out the results of the 

4-parameter functional fit 

(A,t ,6,^) under the assumption 
o 

of zero offsets. Also prints out input 

data values of fitted function and the deviations for each input 
data point. 


Output 

1 
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Output 
# 2 


Prints out 12-parameter fit results, 
treating the 8 S&H offsets as parameters 
to be determined. Also prints out 
Input data, fitted function, and deviations. 


Prints out 4-parameter fit after the 8 
offsets determined above have been subtracted 
from the 32 input data. The 8 offsets 
are printed out as C(3) ,C(4) , . . . ,C(10) . 

Also prints out the offset-corrected input 
data, fitted function, and deviations. Use 
of the line-printer-plotting feature is also useful 
In this 3rd output to provide a 
quick-look check of the results. 

The entire subject of best-fitting experimental data to an expected 
functional form is very complicated and often one uses least-squares or 
related procedures not because they are optimum but because they are readily 
available. Our own use of the function-fitting routine X2MIN(as just described) 
is in this spirit; if the first results of this approach were encouraging, we 
intended later work to verify its suitability. 

The initial results of applying the linefitting offset estimation tech- 
nique appeared very promising. Figure 4 presents the results for SL-2, 

Pass 9, Mode 5, Submode 2, and it is readily seen that the 12 -parameter fit 
is better than the 4-parameter fit assuming zero offsets. Another example 
is provided by Figure 5 showing the resiilts for SL-3, Pass 7/18, Mode 5, 

Submode 2. Another case analyzed was SL-2, Pass 6, M5, SM2, and results from 
these three passes just mentioned are entered in the first three lines of 
Table Z • This table also indicates the averages of the offsets determined 
from these three passes and we were gratified that the pattern of the offsets 
was very similar in these cases. These three cases happened to be among the 
earliest ones we used, and the approach seemed very promising. 

However, there were difficulties soon encountered. Specifically, the 
case of SL-3, Pass 28/39 stopped our apparent progress on the offset-and- 
llnefltting approaches. There was high interest in this particular pass as 
It appeared to be the first pass in the Skylab missions in which the pulse 
compression submodc seemed to be operating and the obvious question was how 
the pulse compression, SMI, and the short pulse, SM2, submodes of Mode 5 compared 


Output 
# 3 
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Figure A. Llneflttlng Results For SL-2, Pass 9, Mode 5, Subiuode 2. 

--Upper Curve Is 4-Parameter Fit With Zero Offsets, Lower Curve 
is 12-Parameter Fit After Subtracting Offsets. 
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Table 2. Results of Linefitting. [Offsets 0^^ In WLlUvolts. 6^^ Defined as “ 0^^ - 0^^^] 
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(PMt>Kl«tioa 

(SL>2) 6 
($1^2) 9 
(SL-3)7/IS 

(SL-3) 28/39 
(SI/>3) 28/39 
(SL-4)79/2* 


Table 2 (continued) « 


ModcCo) 8H 

5 * 

5 2 

5 2 

S 1 

5 2 

5(2) I 


Results of Llneflttlng* 


f fra«tf> 


k 


45 

15.9 

98.9 

12.40 

43 

11.1 

127.4 

422.61 

47 

29.26 

148.4 

428.70 

95 

18.48 

95.56 

>11.61 

47 

24.96 

43.70 

41.94 

94 

25.34 

63.0 

429.0 


e % 

8,38 0.761 

9,19 0.609 

11.07 0.468 

15.07 0.871 

9.49 0.929 

13.27 0.225 

vO 

I 


[A in Millivolts, t^ and 0 in *Nanoseconds , % in Degrees) 
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^gure 6 presents the results of the linefitting with and without offset 
fitting and subtraction for SL-3, Pass 28/39, Mode 5, Subraode 2, and the 
set of estimated offsets is also listed in Table 2, Notice that the average 
of these offsets is more than 10 times as large as in the first three cases 
discussed above. Notice also the considerable difference in risetimes for 
the fits with and without offset corrections for M5, S^]2 of Pass 28/39. 

Proceeding to the pulse compression and applying the same operations 
to the data of Submode 1 of Pass 28/39, Mode 5, the results of Figure 7 
are obtained, with the offsets again listed in Table 2 . Here the diffi- 
culty is even more extreme; the fitting program prefers large negative offsets! 

problem here seems to be that there simply is not ' a sufficient num- 
ber of data points in the earlier portions of the mean return waveform’s 
leading edge, so that the fitting program cannot distinguish between the 
variables 8 and t^ (and also offsets). Figure 8 may be related to this 
point; it shows the various waveform sensitivities to changes in A, t^, B, 
and C as calculated for the particular values of these as determined for 
St-3. Pass 7/18, Mode 5, SM2 fsee Figure 5 ). The function f , itself will 
have the same shape as the derivative (df^/dA) shown in Figure 8 , and it 
is apparent that the (df^/dt^) and df^/dB) curves have somewhat the same 
behavior in the upper half of the "ramp” portion of f^. In the lower half 
of the ramp these two derivatives are clearly distinguishable in their effect 
Qg they have opposite signs. We thinh that the major difference between the 
results for Mode 5, SM2, for SL-3, Pass 7/18 and Pass 28/39 is that the former 
pass had adequate S&H sampling in the earlier half of the ramp and the latter 
pass did not. We think that this may have been due to a greater pointing 
angle which led (because of the change in mean return waveshape as a result 
of Increasing O to a change in tracker bias point such that the tracker posi- 
tioned the S&H gates too late in time, but we can’t prove this conclusively 
from the data of Pass 28/39. 

One other point that might be remarked on Figure 8 is that the sensiti- 
vity to 5 Increases as one moves later in time. This is simply the familar 
result that the later plateau regions are more sensitive to pointing angle 
than is the ramp region where the pointing angle is less than the half- 
bcarawldth of the antenna (when the pointing angle is greater than the half- 
beamwidth, the mean waveform is grossly distorted and no longer characterizable 
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Figure 6. Linefitting Results For SL-3, Pass 28/39, Mode 5, Submode 2. 

Upper Curve Is 4-Parameter Fit With Zero Offsets, Lower Curve 
Is 12-Parameter Fit After Subtracting Offsets. 
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Figure 7. Linefitting Results For SL-3, Pass 28/39, Mode 5, Submode 1. 
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Positlons o£ S&H Gates, H5, SM2. 



Figure 8 . Theoretical waveform derivatives calculated for A=148.4rav, t^-28.7ns, 
g«ll. 07ns, and degrees, as found by 12-parametcr fit to 

Pass 7/18 of SL-3, Mode 5, Submode 2. S&H gate positions in each 
subsubmode are indicated by arrows. 
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aa having a ramp and a plateau). Thus for determining the pointing angle 
from the S&H waveform data (as needed in O work for example) , the 100ns 
pulsewidth submode would seem preferable to the 20ns submode* 

Figure 9 shows again the S&H offset estimation via linefitting for 
SL-4, Pass 79/24, MS, SMI. Here a different problem is built into the again 
unacceptable results. It is now thought that the antenna cliaracteristics 
changed between SL-3 and SL-4 and without adequate knowledge of the new 
antenna patterns, the fitting routine’s cpnstant y is not adequately known. 
[Recall that y is a beamwidth-related factor determined as specified by the 
earlier Applied Science Associates Report on Task D as cited in the discus- 
sion of f(t) several pages earlier.] In addition to the uncertainty in y, 
the results of Figure 9, with offsets as listed in Table 2 , have very large 
positive offsets. These offsets show an apparently quite different distri- 
bution about the mean offset than all other entries in Table 2 [see specifi- 
cally the 6^ in that Table]. 

Finally, there is nothing in the preceeding discussion which restricts 

wv •Duuuiuutiia j- “v/j- ^ wiaw - 

to 4 subsubmodes and to 32 data points, and the general linefitting procedure 
has also been applied to three different submode 0 results from different M5, 
SL-4 passes. These again are subject to the uncertainty in the beamwidth- 
p^iated constant y. The results from these passes are presented in Figures 
10, 11, and 12, and here again the results are only partly consistent. There 
is very little more that can be said about the S&H offset problem at this 
time; a general approach has been described and the results are mixed. Some 
jiioxe variations on this should be carried out, and we again emphasize the 
view that the present report is only a statement of progress to date on an 
unfinished problem. The remaining pages in this section are addressed to a 
different waveform-related problem, tliat of the mean return waveform’s statis 
tical properties. 
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Solid line = fitted curves 



Figure 11. Linefitting Results For SL-4, Pass 83/29, Mode 5, Submode 0. 

Upper Curve Is 4-Parameter Fit With Zero Offsets, Lower Curve 
Is 12-Parametcr Fit After Subtracting Offsets. ' 
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Figure 12. Linefitting Results For SL-4, Pass 86/32, Mode 5, Submode 0. 

Upper Curve Is 4-Parameter Fit With Zero Offsets, Lower Curve 
Is 12-Parameter Fit After Subtracting Offsets. 
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Most of the discussion thus far has been concerned with the concept 
of the mean return waveform, and the arithmetic averages of S&H gate reading 
have been used as estimates of this mean return waveform- It is useful 
to perform the sum of squares of gate readings at the same time standard 
In order to be able to estimate the standard deviation as well as the mean, 
and we practically always do provide standard deviation estimates in all 
our S&H-r elated programming. 

But the mean and standard deviation are not sufficient to characterize 
an unknown distribution, and one needs additional information. In the 
following discussion we will Ignore any tracker time-jitter effects by 
assuming that either such effects are negligible or have already been cor- 
rected for (by the time-realignment process for example) before the histo- 
gramming procedure is applied. 

By selecting several points (at different times) on the mean waveform 
and producing a histogram at each point of the frequency of occurrence of 
"each '"voirage , wc can gcc -ac eacn-oz 'ciiese-pciivcs a coarse estxmate oi. 
the probability distribution function (pdf) of the S&H voltage. The ”coarse- 
ness" of the pdf estimate is reduced as one increases the number of histo- 
gram bins but this cannot be done without limit because of the requirement 
of having an adequate number of samples in each bin- For this reason and 
because of practical limitations on how much computer space is to be tied 
up or how many output points have to be plotted, an upper limit of the order 
of 100 bins is reasonable for at least the first try at the problem of obtain 
Ing the waveform statistical properties via the histogram approach. 

A general-purpose histogramming subroutine HIST has been written for 
this purpose.* An entry point HISTI is provided for initialization and 
another entry point HISTO prints out the results of histogramming. We will 
sketch here the use of HIST, HISTI, and HISTO for the Individual Inputs to 


*At the time of this report's writing, this histogramming procedure had not 
yet been implemented. The subroutines HIST had been debugged and delivered 
to Wallops Flight Center and the first results were being obtained while 
this report was being prepared. 
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S&H gates 1,3,5 and 7 In each subsubmode of each submode of Interest* For 
a first trial, the histogramming range should be from -0.050 volts to 40.450 
volts, virtually the entire non-saturated range of the S-193 S&H gates. 

We present immediately below the outline of the form a main program will 
have to carry this out, under the assumption that all S&H data in each given 
subsubmode will be used to achieve the maximum precision available for esti- 
mates of the statistical properties at S&H gates 1,3,5, and 7 on the mean 
return waveform. Appendix C provides a program description for subroutine 
HIST: it will be necessary to refer to that description of HIST to under- 
stand its use in the "main program" below, 

**Main Program” 


other DIMENSION stuff _ 

DIMENSION JHl(lOO) ,JH2(100) , JH3(100),JH4(100) ,SUMSl(2) ,.SUMS2(2) 

DIMENSION SUMS3(2),SUMS4(2) 

DATA NB/100/,NSIGMA/1/,XBAR/. 250/, SIGMA/. 200/ 

Other initialization 
^ ^ new submode 
I — i»>.new subsubmode 

I j initialization within suDsubmoue 

I CALL HISTI(XBAR,SIGMA,NSIG,NB,WB,XL,JH1,NHE1,SUMS1) 

I I CALL HISTI(XBAR,SIGMA,NSIG,ITB,WB,XL,JH2,NHE2,SUMS2) 

I I CALL HISTI(XBAR,SIGI4A,NSIG,NB,WB,XL,JH3,NHE3;SUMS3) 

* * CALL HISTI(XBAR,SIGMA,NSIG,NB,V7B,XL,JH4,NHE4,SUMS4) « . 

I j j ^oop to read S&H data, 104 sets per frame, for all S'^ frAifteS^ ' 

, I read S&H //I into SHI 

I I read S&H //3 into SH2 

I I read S&H #5 into SH3 

I read S&H #7 into SH4 
1 I j CALL HIST(SH1,WB,XL,JH1,NHE1,SUMS1) 

I I I CALL HIST(SH2,WB,XL,JH2,NHE2,SUtIS2) 

» » CALL HIST(SH3,WB,XL,JH3,NHE3,SUMS3) 
j I J CALL HIST (SH4,WB, XL, JH4,NHE4, SUMS ) 

I L*. -.-end of loop for individual set of S&H values 
I 1 write heading for S&H //I 
I CALL HISTO(NB,WB,JH ,NHE .SUMS ) 

j write heading for S&H //2 

* CALL IIIST0(NB,WB,JH2,NHE2,SUMS2) 

I I write heading for S&H //3 

I I CALL HIST0(NB,WB,JH3,NHE3,SUMS3) 

* • write heading for S&H M 

I j CALL HIST0(NB,WB,^4,NHE4,SUMS4) 

I — go to next S M 

I go to next SM 

finish out problem, etc. 
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SECTION III. WAVEFOKM TIt4E- REALIGNMENT AND AVERAGING 

The background for the results of this section is provided by an earlier 
Applied Science Associated report ["Final Report on Task A - Engineering 
Studies Related to the Skylab Program", prepared under Contract NAS6-2307, 
approx. July 1973]; that report is also the source for the weighting func- 
tion used to estimate true altitude tracker position, given the average of 
the tracker’s positions over an approximately 1/8 second time interval. We 
%dll first discuss briefly the need for time-realignment and the procedures 
to carry out a realignment followed by an averaging operation, then repre- 
sentative results will be presented. The computer program details involved 
in this are in Appendix D. 

Because the 8 S&H gates are fed directly to the digital delay generator 
(DDG) in the S-193 radar altimeter's altitude tracking loop, finite altitude 
rates and/or tracker jitter will lead to each S&H set’s sampling different 
time-positions on a return waveform. Simply averaging together aU measure- 

Bents by S&H gate #1, gate #2 8 to form average sample points 91 . 92 . 

without correcting for the tracker jitter shouid ieaa ro a mean wave- 
form distortion and an increase in the variance of the 8 sample points. A 
more detailed investigation of jitter effects on waveform and variance is 
provided by another Applied Science Associates Report ["Final Report on 
Task D - Engineering Studies Related to the GEOS-C Radar Altimeter by 
L. S. Miller and G. S. Brown, prepared under Contract NAS6-2307, May, 1974]. 
Briefly, the mean waveform is "smeared”, essentially by a convolution with 
the tracker jitter process’s probability density function and there is an 
attendant increase in variance estimates. The waveshape effects should be 
negligible for the S-193 100 nanosecond pulsewldth modes but significant 
for the short pulse submodes of Mode 5, the 10 ns (nominal design pulsewidth) 
direct or pulse compression submodes. 

The correction for the tracker jitter involves comparing Instantaneous 
tracker position with "true" tracker position; this difference is a measure 
of how far the 8 S&H gate results must be repositioned earlier or later in 
time before adding up a number of separate sets of S&H data to determine an 
average waveform over a suitable averaging period. There are two unknown 
“quiitlties to be estimated to find this instantaneous tracker error however. 
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the "true tracker position" (or equivalently, the "true" altitude) and the 
Instantaneous tracker position. 

To determine the "true altitude" as a function of time, a low-order poly- 
nomial is least-squares fitted to the altitude data over an entire submode - 
This is in effect a type of smoothing operation. In practice, the maximum 
degree is 4 for the polynomial fitting [done in subroutine POLRG as explained 

In Appendix D]. 

The altitude outputs from the S-193 altimeter at the 8 per frame rate 
are not 8/frame measurements of the tracker's instantaneous position. Rather, 
each altitude output is an average of the tracker's position during the ~l/8 
frame preceedlng that output. Part of the Applied Science Associates Task A 
Report, July 1973, was devoted to deriving a weighting function to estimate 
tracker instantaneous positions, and that weighting function is -used in this 


work. 


The 8/frame instantaneous tracker position estimates are used to produce 
the needed 104/frame S&H time corrections, and a spline is used for this 


[sC' 


t . .t j T»1 _ 


time— realigned sets of Cull readings arc sui 


a set of time bins [as described in the Applied Science Task A Report and 
summarized by the time sketch in subroutine VTADD in Appendix D of this 
report], and the results are written on tape on a frame-by— frame basis. A 


subsequent averaging operation reads the frame-by-frame tape, and prints and 
plots the results on a 10 frame basis [this last program was developed by 
Wallops personnel and is not described in this report]. 


Figure 13 summarizes the overall time-realignment and averaging as just 
described. Figure 14 shows a portion of actual altitude data as it comes out 
of box D of Figure 13 ; also indicated are the 8/frame instantaneous tracker 
error estimates out of box E of Figure 13. Then Figure 15 repeats these 8/frame 
Instantaneous error estimates and shows the 104/frame individual tracker 
errors as produced by the spline - these 104/frame error estimates provide 
the time-corrections to be applied before summing the 104/frame sets of S&H 

data. 


It is necessary that the S&H data already be corrected for dc offsets 
before carrying out the averaging after the time-realignment. This is because 
each time bin's contents includes contributions from several different S&H 
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S-193 SL-2, Pass 9, Mode 5, Submode 1, Subsubmode 0. 
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gates over the time period of a frame, and there is no way to calculate or 
correct for this combination of several different S&H offsets after the 
time-realignment is applied. 

(Me additional problem (not shown on Figure 13) is that it is possible 
that some of the S&H readings in the primary data tape (from JSC in Houston) 
may fall outside the allowed region of the S&H calibration curves; this occurs 
because of the way in which the calibration curves were handled in the JSC 
processing. To avoid contamination of the waveforms by these invalid data, 
a test procedure has been built into the programs at Wallops to detect the 
occurrence of S&H readings outside the allowed range, and to flag these data 
so that they are not included in the waveform averages determined in the 
later program steps. Consequently there will not necessarily be 104 contri- 
butions per frame to the waveform averages, but 104 will be the maximum num- 
ber possible per frame. 

Decisions have yet to be made about how to handle date for which there 
are appreciable numbers of loss-of-lock indications in the altitude tracker ^s 

^ ^ j gwi r»-F dalhj* AT»d « oS 

■"OUC'pUi.." - it '"Jio '•t'o-'-*-'*- "'—“—'to'*'* ' >-—• — . 

already mentioned in the introduction to this report, no suitable set of 
data quality criteria has yet been found which would permit the implementa- 
tion of an automatic data editing part of the overall waveform processing. 

Figure 16 presents a preliminary indication of the time-realignment re- 
sults from SL-2 Pass 9, Mode 5, submode 2; what has not been shown on the 
figure is that points shown are of unequal weights. The next step would be 
to form appropriate weighted means and then use these as input to a final 
pass through X2MIN to determine final mean return waveshape parameters. It 
is apparent however that no striking change in the leading edge risetime has 
appeared - and this is the parameter of the greatest Interest in attempts to 
detect sea-state from S-193 waveform measurements. Another disturbing feature 
not displayed in the figure is that no changes were found in variances esti- 
mated from the time-realignment procedure as opposed to variances from simply 
summing up all of the measurements of a specified S&H gate - we had expected 
to see variance changes approaching 25%. 

Figure 17 presents results for the time-realignment process applied to 
SL-3 Pass 28/39, Mode 5, Submode 1; no offset corrections have been applied 
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Figure 16. SL-2, Pass 9, Mode 5, Sufcmode 2 Time Re-alignmerit Results 



No offset corrections have been applied. Results shown are averages over 10 frames 
In each sub-submode (except sub-submode 2: for which only 5 frames of data are taken). 
Note that the Individual points below fcr the time-realigned waveform are not of 



Figure 17. SL-3, Pass 28/39, Mode 5, Submode 1 (Pulse Compression) Time-Realignment Results. 
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to the S&H data before time realignment. This figure displays one of the 
hoped-for outcomes of the time-realignment as there seem to be data which 
are recoverable in the leading edge earlier than the first S&H gate. Since 
the S&H set started too late on the waveform, even in subsubmode 0, these 
data are very deslreable. However, we again here have totally unknown S&H 
offsets. The only chance of recovering the desired information might be 
some sort of iterative approach using first the offset-determining procedure, 
then the time-realignment; by going back and forth, and carrying lineshape 
parameters from the time- realignment back to the offset determination, it 
might be possible to recover something. We have not yet tried this. 

Returning to the SL-2 Pass 9, M5, SM2 time-realignment results, no sig- 
nificant changes were obtained for a variety of different guesses at S&H 
offsets. Eventually growing frustrated by not being able to make the situa- 
tion better, we tried to make things worse just to determine that there were 
no overlooked sign errors. For instance, first we replaced the weighting 
function by a single unit weight (which is equivalent to simply bypassing the 
box E of Figure 13) - again, no significant changes were found in the results. 

ly we reversed the sign ot the time— error correction wirh che results 
shown in Figure 18. There it does appear that the risetime is shorter for 
the non-time-realigned results than for the deliberately erroreously corrected 
results . 

These results all indicate that we are not now gaining enough informa- 
tion from the time-realignment to justify applying it in any automatic or 
routine fashion to the remainder of the S— 193 waveform data. Some additional 
small-scale investigation of time-realignment should be continued on a case- 
by-case basis, and the work presented in this report provides the basic build- 
ing blocks for this. The change in the short-pulse mode's pulsewidth from 
nominal design width of 10 ns to an actual pulsewidth closer to 20 ns is 
probably the largest single factor contributing to our failure to realize 
benefits from time-realignment. It is worth noting that our S-193 experience 
further supports the conclusions of the earlier Applied Science Associates 
Report (our Task D of Contract NAS6-2307] which indicated that waveform tlme- 
reallgdment was not going to be necessary for GEOS-C if that system met its 
specifications. 
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APPENDIX A. SUBROUTINE X2MIN FOR GENERAL FUNCTION FITTING TO DATA. 


In our investigation of the S&H offset in waveform data from the S-193 
radar altimeter, we used a general-purpose subroutine X2MIN for fitting a 
theoretical function to the experimental data. X2MIN is sufficiently general 
as to merit further description in this appendix. We did not Inveiit X2MIN, 
and the version presented here differs only in very minor details from a 
program obtained from Professor Dwight W. Carpenter at the Duke University 
Physics Department (Dr. Carpenter is now at Limestone College, Gaffney, 

South Carolina). Dr. Carpenter’s program was In turn an extensive revision 
of an earlier general-purpose fitting routine from the University of Illinois. 
The following description is entirely derived from Dr. Carpenter’s program 
write-up at Duke University. 

X2MIN is a FORTRAN subroutine to fit an arbitrary function to given 
data points by minimizing the weighted squared deviations of the points from 
the function. The function is of the form. 


FX(y^,X^,.. 


V A A 






The X’s are the independent variables, the A’s are the variable parameters 
to be adjusted by the program, and the C’s are constant parameters. The 
data consists of a set of NP points of measured at given values of 
associated error estimates cJY^; i = 1, NP. 

„ T [Y^ - FX(X^, A.C)]^ 

1=1 

The program requires initial guesses for the variable parameters A. Succes- 
sive corrections are made to the A’s until a relative minimum in x is reached. 
Reasonably good initial guesses are often necessary. Additional constraints 
may be placed on the A’s by providing an error estimate oA on the Initial 
value of A^. For each constrained A^, the term 

t<4itial-4inal>'/<«*">'’ 

2 

is added to the x • 
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The Iteration continues until the fractional change in X step 

is less than a given limit (CHILIM) , or until a given maximum number of 
interations (MAXITR) is reached. 

The user must supply a function type FORTRAN subroutine^ FUNCTION FX, 
which will evaluate FX (and, optionally, its derivatives with respect to the 
A*s) for a given set of X*s, A*s, and C*s. The data deck organization for 
any one problem is as follows: 


Card Type 

Information 

No- 

Short Cuts 

Comment 

Any 

1 

Repeat? 

Size 

NP ,NX,NA,NC , CHILIM,>IAXITR,NPLOT 

1 

+ 

Data Label 

X,Y,OY. names 

1 

4’ 

Repeat? 

Data Point 

X,Y,(3Y values 

NP 

4 

Var. Param (A) 

A name. Initial value, step, constraint 

NA 

4 


Const Param (C) C name, value 



This may be followed by succeeding problems.* The Short Cuts may be used on 
succeeding problems if portions of the input data are to be repeated. A 
single card ^ich preceeds the first of the problem decks sets up a set of 
ten general flags INFIAG(IO) which are available through the labelled COMMON 

area, COMl. 

The calling procedure for X2MXN is as follows : 


— Misc Job Control Cards — (Input on device 1, output on 3) 

COMMON /COMl/ INFLAG (10) 

READd, 5) (INFLAG(I) ,1=1,10) 

5 FORMAT (1015) 

CALL X2MIN 

STOP 

END 

Subroutine X2MIN calls the additional subroutines LSQMIN, ^^INV20, 

XPROB, and PL0T4 (with entry points PLOTWD, PLOTS, and PL0T3L) as well as 
function FX. Following details of the X2MIN input data and of FX in general, 
source listings are provided for X2MIN, LSQMIN, MINV20 and XPROB. In addi- 
tion to facilitate use of X2MIN elsewhere and to ease debugging or modifi- 
cation, logical flow diagrams are provided for X2MIN and LSQMIN. Because 
the details are computer— dependent, no source is provided for PL0T4 (with 
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PLOTWD> PL0T3, and PL0T3L) ; however a description of PL0T4 is provided. 
X2MIH Input Data 

Six types of data cards are described in the order of their appearance 
In the data deck: 

1. COMMENT CARDS 

The comment cards will si^Tply be read and printed out. At least one 
comment card must appear. Eat card except the last must have .... (four 
periods) punched in columns 77-80; the last card must not have them. 

If any comment card startS^with the word REPEAT in columns l“6f the 
program will expect the Variable Parameter (A) cards to .follow the last 
comment card. The Size, Data Label, and Data Point information will be 
carried over from the preceding problem. 


2. SIZE CARD. _ 

Seven numbers appear on ®ne card in successive fields of 10: 


NUMBER 

FORMAT 

DESCRIPTION 

ALLOWED VALUES 

UP 

5X, 

15 

Number of data points 

1-100 

NX 

5X, 

15 

Number of independent variables (X’s) 

1-10 

NA 

5X, 

15 

Number of variable parameters (A*s) 

1-20 

NC 

5X, 

15 

Number of constant parameters (C*s) 

0-20 

CHILIM 

5X, 

F5.0 

2 2 

Ax /X convergence test. (If left 
blank, 0.001 is used) 

Any positive value 

MAXITR 

5X, 

15 

Maximum ^number of iter act ions taken. 
(If lef£- blank, 30 is used) 

Any positive integer 

NPLOT 

5X, 

15 

1 causes plots to be deleted 

0 (blank) -1 


3. DATA LABEL CARD 

Alphabetic characters, 10 columns each, giving the names of the variables 
Y, aY (NX + 2 names) in sequence. These will be used as column 

headings in the printout. 


If the word REPEAT appears in the first 6 columns (instead of the labels) , 
the Data Point Cards must bcj' deleted. The Data Label and Data Point informa- 
tion will be carried over from the previous problem. 
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If the words STAT ERROR are inserted on the Data Label Card as the 
label for 0Y, the value taken for oY Is (not /y” ) for each point. A 
lower limit of 0.1 is arbitrarily Imposed on FX for this calculation. The 
OY columns on the Data Point Cards are ignored in this option. 

4. DATA POINT CARDS 

One card for each data point as follows: 

X^, ^ format (8KL0.0). If the OY are left blank (and 

the STAT ERROR is not specified) , all points will be given equal weight and 
portions of the output dealing with errors will be deleted. 


5, VARIABLE PARAMETER CARDS (A*S) 

One card for each parameter as follows: 


Parameter Name 

(AlO) 

Alphabetic name of parameter 

Continuation Flag 

(Al) 

If C, the Initial value is taken as the 
final value in the preceding problem. 
If blank, the new initial value is 



taken -fmm the .nerd r 

Initial Value 

(F14.0) 

Initial Value 

Step Size 

(F15.0) 

Step (AA) used to find dY/dA, or blank if 
subroutine computes dY/dA directly. 

Constraint 

(F15.0) 

Standard deviation (OA) assi^ed to Initial 
value. If blank (or zero), no con- 
straint is made. 


6. CONSTANT PARAMETER CARDS (C*s) 

One card for each (if any) constant parameter as follows . 

Parameter Name (AlO) 

Continuation Flag (Al) — as above in A*s 
Constant Value (F14.0) 

evaluation of derivatives 

The program must evaluate the derivatives of FX with respect to the A 
parameters. This may be done in one of two ways for each parameter: 

(1) A step size AA (non-zero) may be given on the A cards. The program 
then makes a step AA in A and takes AFX/AA as the derivative. The step size 
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nust be chosen small enough that the derivative Is essentially constant over 
the step, but large enough to make a significant chainge in FX, An appropri- 
ate step size is often critical. It may be checked by varying the step size 
to see that the results are independent of reasonable changes. 

(2) A formula 'which evaluates dFX/dA for each A may be wltten into 
the function routine FX. In this case, the function must return 

FX » dFX/dA(JA) whenever JA > 0 

and 

FX « FX function whenever JA = 0. 

JA is a flag given the subroutine by the main program. The step size field 
must be left blank (or zero) for each parameter whose derivative is to be 
computed in this way. 

The two methods may be intermixed if desired. The main program gives 
the flag JA > 0 for those A not given step sizes. 

PRINTOUT 

The printout contains the following Information for each problem: 

All input information. 

2 

The value of x the A’s at each iteration. 

/ ’ 7 — 

The final errors on the A’s, /<(AA) > 

2 

The degrees of freedom and X“ probability (high number for good fit, 
low for bad, in percent) 

The value of the best fit function FX at every data point 

The deviation of each data point from the best fit FX 

2 2 
The (devlatlon/aY) for each point, i-e., the contribution to x 

The number of the worst point 

A plot of the deviations in Y^ vs X^ 

A plot of the Y^ vs X^ 

A correlation matrix 

The error matrix is (EM)^^ *» <AA^ • AA^> where AA^ is the error in the 
parameter A^. The rms error on A^ is /<AA^ • AA~> . The correlation matrix 
is related to the error matrix by: 

(CM)^^ » (EM)^j/ /:AA^>. /<AAj> . 
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(CM)^j will range from 0.0 for uncorrelated parameters to + 1.0 for com- 
pletely correlated parameters. 


ERROR MESSAGES 

The message SINGULAR MATRIX will occur if a matrix is encountered which 
cannot be inverted as required. This usually indicates an error in the data 
or in FX. 


FUNCTION FX 

The function subroutine FX must be supplied. The FUNCTION, COMMON, and 

DIMENSION cards must be as shoxra. The routine must evaluate FX for point JP; 

FUNCTION FX (JP,JA) 

COMMON/COMI/INFIAG (10) 

COMMON X,Y,A,NA,NX,NC,C,NP 
DIMENSION X(100, 10) ,Y(100) ,A(20) ,C(20) 

——test on JP^O 

FX ^ function evaluated for X(J,1),. ..X(JP,NX) ,A(1) . .A(NA) , 
C(l),.. .C(NC) 

, ' RETURN 
END 

In X(JP, JX), JP is the point nunber (l£JP<OT) and JX is the independent 
variable number (1^JX<NX) . Immediately upon entering X2MIN for the first 
time, there is a statement *'F = FX(-1,~1)'*; after the first problem is finished, 
the program will execute the statement ”F FX (0,0)” i mme diately prior to 
reading each new problem. The values -1 and 0 for JP are provided for possible 
Initialization in FX, and must be tested for. The use of COMl is optional; 
only if some use is to be made of flags INFLAG(IO) must COMI be provided in 


If the function is to evaluate derivatives explicitly (as described above 
as method 2) JA must be tested; 

IF (JA) 1, 1, 2 

1 FX “ FUNCTION 
RETURN 

2 FX * derivative dFX/dA(JA) 

RETURN 

END 

example 


P is measured for various T and V for gasses with known R. One wants 
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to find AA and BB in the following equation: 
P = R*T/(V - BB) - AA/V^ 

The FX subroutine might be: 


FUNCTION FX (JP,JA) 

COMMON X,Y,A,NA,NX,NC,C 

DIMENSION X(100,10) ,Y(100) ,A(20) ,C(20) 

IF (JP.LE.O) RETURN 
T=X(JP,1) 

V=X(JP,2) 

R«»C(1) 

AA=A(1) 

BB«A(2) 

IF(JA) 1,1,2 

1 PX-R*T/(V-BB)-AA/V**2 

return 

2 GO TO (3,4) ,JA 

3 FX=-1.0/V**2 
RETURN 

4 FX=R*T/(V-BB)**2 
RETURN 

END 

A typical data deck for tnis problem mighc then appear lollows (notice 
the use of the REPEAT feature, and the continuation flag C on parameter BB 
in the repeated part) : 


■5 s 

? 7 

0 

CM 

1 

H 25 

*- 30 

0 

•it 

1 

f T 

FIT AA AND 

t 

BB FOR GAS 

T 1 

1 

NP=* 5 

NX= 2 

NA= 2 

NC= 1 

TEMP 

VOL 

PRESS 

P ERROR 

300.0 

27.6 

18.4 

0.10 

230.0 

20.0 

32.4 

0.15 

240.0 

20.0 

15.2 

0.10 


AA .03 

BB .01 

R 50.2 


REPEAT WITH AA CONSTRAINED TO .03 4-. 01 
AND USING STEP SIZES TO FIND DERIVATIVES (CONTINUE BB) 
AA .03 .0001 .01 

_ BBC .001 

R 


n 


o 

CO 




50.2 
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NOTES 

The constraint in the A’s is useful for two things: 

(1) Introducing outside knowledge about a parameter 

(2) keeping a variable from running wild in the early stages of 
fitting* The constraint may be released in later stages - 

The function FX has available to it NA and NC which may be useful in 
having the number of terms set at run time, or changed during a run. 

The C’s may be used as program switches by the function FX as well as 
actual numerical constants. Thus several different functions could be In- 
cluded in one FX. The flags INFLAG(IO) are provided for additional switches, 
but the first three of these are used by PL0T4 as described later. 

If the step sizes are used, one should see what the effect of the choice 
of size has on the results. 

A high correlation coefficient between two parameters (|c|>-9) means 
variations in the two parameters have a similar effect on the function, 

j.M. ^ ~ T . ^ .... ^-4 .. .--I ^ .-"i. /-. 

'ItajLU. tw ...... -w.vk^ 

should try to find a parameterization which will give low correleations. 
Beware uf C>.98. 

The following material will be presented in the order listed here: 

X2MIN source listing, X2MIH flow diagram, LSQMXN source, LSQMIN flow diagram, 
MINV20 source listing, XPkOB source listing, and description (only) of 
PL0T4. 
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00 113 

113 CIUJkM,*) * 


If 00 TO 101 
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f K -«i ■ — 
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nu- n -a — — - — — 
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FL0T4, the Line-Pointer-Plot ting Subroutine Called by X2HIN ' ' 

As now used, X2MIN calls a line-printer-plotting subroutine PL0T4 
which has additional entry points PLOTWD, PL0T3 and PL0T3L.. Since PL0i4 
foms a plot image by storing one character in each byte (8 bits) of a 
4-byte word, the subroutine is specific to the IBM 360/370 series com- 
p^iter we used^ for this reason we do not provide here the source listing 
for PL0T4, However, a brief description follows for those users who would 
wish to Incorporate line printer plots at their own computer • (The alter- 
native to this would be to delete all statements in X2MIN which call PL0T4, 
etc.; that is, delete all statements between 209 and 579 in the X2MIN listing. 
Statement 209 must be left in, 579 must be removed, and the ”go to 6" state- 
ment Immediately following 579 will then appear immediately after statement 
209 to provide the necessary transfer of control.) 

The size of the plot is set up automatically at the first call of PLOTS, 
additional calls to PL0T3 or PL0T3L provide additional data to the plot,, and 
the call to PL0T4 prints out the accumulated plot image and resets it. 

The automatic setup provides that all the points from the first PLOTS call 
are included within the axis limits. The automatic setup sets 10 print 
spaces per grid line in the x-dimension and 5 print spaces per grid line in 
the y— dimension with the grid lines adjusted to give round or nice numbers 
and hence the points plotted won’t in general extend to the edges of the plot 
space. The maximum grid widths allowed for NGX (the number of X— grid lines) 
and NGY (Y-grid) are both preset to 12, as appropriate to a line printer, 
but a call to PLOTWD (MAXGX,MAXGY) will change these limits to MAXGX (maximum 
NGX) and MAXGY (maximum NGY). Reducing MAXGX to 7 permits output plotting on 
the 80-character line of Teletype-like devices. 

In X2MIN, the labelled common area COMl provides ten flags, INFLAG ( 10 ) , 
for general setup; the first three of these are used to set the type of out- . 
put printer plots and the plot size. INFLAG (1) _> 3 produces printer plots 
of both of the fitted function with experimental points superimposed and of 
the deviations between experimental and fitted points, INFLAG (1) = 2 provides 
only the plot of the deviations, INFLAG(l) = 1 plots only the fitted function 
with experimental data, and INFLAG(l) ^ 0 suppresses all plots. If XNFLAG(2) 
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0, INFLAG (2) and INFLAG(3) carry new values for MAXGX and MAXCY in a 
PLOTWD(MAXGX,MAXGY) call; otherwise MAXGX and MAXGY are both 12. 

FL0T3(C,X,Y,N) plots the character C at the N points X,Y(X and Y are 
each one'- dimensional arrays) ; the plot character C replaces whatever was 
previously in the plot image. 

FL0T3L(C,XX,YY,N) is called with either 0 or 1; if N=l, a line is 
drawn (using character C) to the single point XX,YY from the point input of 
the previous PL0T3L call, provided that N=1 in the previous call also. The 
first PL0T3L call with N-1 provides no output, and a N=0 call cuts the string 
to permit a new "first” N=1 call. 

PLOT4 (LX, XLAB,LY,YLAB,LT, TITLE) prints out the accumulated plot image 
and resets it. XLAB, YLAB, and TITLE are character strings for the x and y 
axes and the title, with LX,LY, and LT being the string lengths. 
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APPENDIX B. SUBROUTINE FX. 

The following pages provide the source listing and the flow diagram 
for subroutine FX which supplies the llneshape function and its derivatives 
with respect to each of the 12 lineshape parameters; the mathematical de- 
scription of this FX has already been given in Section II. FX is for use in, 
and follows the rules established by, the general line-fitting subroutine 
X2MIN as described in Appendix A. 

Because the Bessel Functions, and I^, and the probability function 
arc needed in FX, listings are provided for the Bessel function subroutine 
BESI and the probability function subroutine NDTR. Both BESI and NDTR are 
taken directly from the IBM Scientific Subroutine Package and are described 
more fully there. 



pyvcTro*! FXMP,JM 
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appendix C* SUBROUTINE HIST 

The following several pages provide the description of the subroutine 
HIST (with entry points HISTI and HISTO) which is intended for use in 
obtaining the probability density function at one or several points on the 
mean return waveform* A brief description of how to use HIST for this has 
already been provided at the end of Section II of this report. 



- 66 - 


Subroutlne Name and Argument List 

Subroutine HISTI (XBAR,S IGMA,NSIG,NB ,WB , XL, JH,NHE , SUMS) 

Entry HIST(XI,WB,XL, JH,NHE,SUMS) 

Entry HlSTO(NB,WB,WL,JH,NHE,SUMS) 

Description 

Subroutine HISTI, with entry points HIST, HISTO, is to allow the histo- 
gramming of a number of observations of the variable XI. A call to HISTI 
initializes the subroutine and calculates the width and location of individ- 
ual histogram bins. Calls to HIST generate the histogram desired, and a 
final call to HISTO causes print-out of the results; the call to HISTO does 
not disturb the histogramming storage locations so that one may display an 
Intermediate result by HISTO and then continue to build up the histogram 
by more calls to HIST. All variables used are listed in calling lists so 
that it is possible to have several histograms being built up at the same 
time* 

Fm present subroutine dimensions, a maximum of 100 histogram storage 
locations is allowed; these dimensions are easily Increased if desired. 


Subroutine Argtjments (and Dimensions) in Calling List 


XBAR \ These are input arguments to establish the range which 

SIGMA I the histogramming will cover uniformly with a total of 

MSIG ( NB histogram bins. XBAR is the value at which the 

NB I histogramming routine is centered, and (NSIG*SIGMA) 

Is the width of the total range to be histogrammed. 

Input data lying outside the range XBAR+(NSIG*SIGMA) 
will be entered either in bin //I or bin //NB, depending 
on whether the input is below or above the range Indicated. 


WB ) These values are calculated in HISTI, given the inputs 

XL j XBAR, SIGMA, NSIG, and NB. WB is the width of an 

individual bln, and XL is used in index computation in 
HIST- 


JH(IOO) Integer array JH is the histogramming "count" storage. 
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NHE 

SUMS (2) 


XI 


Counter recording total number of histogram counts in 
the locations of JH (100) . 

Two double-precision variables used for estimating an 
overall mean and standard deviation for the histogrammed 
variable XI; SUMS(l) is a running subtotal of individual 
input values XI, and SUMS (2) stores the subtotal of the 
squares of XI. 

Individual input value of the variable to be histogrammed. 


Tnput/Output Performed by Subroutine 

No input except through calling list. HISTO writes (on unit 3) the 
probability density function estimates obtained from the histogram, together 
with estimates of the mean and standard deviation. The pd.f. estimates may 
be converted to total histogram-bin counts, if desired, by multiplying all 
iwl.f, estimates by NHE, the total number of histogram entries. 
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Program Listing for HISTI»HIST, and HISTO 


SUBROUTINE HI STI (XBAR^ SI GMA^ NSI VE>XL^ JH>UHE> SUMS ) 
DIMEN S 1 Crj JH ( I 0 0 ) ^ XPR’JT ( 1 13 > # PPR;'J T C 1 0 ) > SUM S ( 2 ) 

DOUBLE PRECISION SUMS/ DXMN# DSG^ DQ 

IF <NB.GT. 100) Nb=l 00 

MHE=0 

SUMS(1)=Q.D0 
SUMS(2)=0-D0 
DO 10 I=1/NB 
10 JHCI)=0 


Q=NB-3 

Q=NSIG»SIGMA/Q 

VB=Q 

XL=XBAR-<NB-2)*Q/2. 

RETURN 

ENTRY ;HIST<XI/ VB/XL/ JH/NHE/SUMS) 
Is2+IMTC <XI-XL)/VB) 

IF (XI.LT.XL) 1=1 
IF (I.GTtNB) I=iIB 
20 MHE=NHE+1 

JH< I )=JH< I )+l 
SUMSC 1 )=SUMSC 1 )+XI 
SUMSC2)=SUMS(2)+XI*XI 
RETURN 

■ENTRY. HIS.TCCNBv XL,,* ^ 

QaNHE 










■is 




DQ=Q 

DXMN= SUMSC 1 ) /DQ 
DSG=SUKS(2)-DG-» dxmn » dxkn 
SG=DSG/(DQ- I .DO) 
SG=SGRTCABSCSG) ) 

XMN= DXMM 
X=XL-WB/2. 


1 = 0 

WRITE (3/30) NHE/VB 

30 FCRMATC/^ FCLLCVING IS H I STOGRaM-BEHI VED PDF FOR / IS/ 

HISTCGRAI^ ENTRIES/'/' IN FORM ( CENTER/ PDF ) . BIN WIDTH 


'/E13. 6 


45 J=0 

DC 35 K=l/4 
1=1 + 1 

IF C 1 .GT.NB) GO TO 50 


J=J+l 

XPRNT(J)=X 

X=X+Wb 


QQ^JHC I > 

35 PPRNTC J ) = QQ/Q 
50 IF CJ.GT.O) WRITE C3/40> 
40 FCRMATC ' '/ 4C ' C '/ F3 . 3/ '/ 

IF (l.LT.NB) GC TO 45 
WRITE (3/55) tJHE/XMN/SG 
55 FCRMATC' FCH ABOVE'/ 15/' 
I E13.6/) 

RETURN 

END 


(XPP.NT(K) /PPRIJT(K)/ X=l/J) 

'/F7. 5/ ') ') ) 

£NTRIES/MEA:J= '>E1 3. 6/ '/ AND STD DEV= 




Jt 
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appendix D. PROGRAMS FOR TIME- REALIGNMENT AND AVERAGING. 

The various program descriptions in Appendix D carry out details of 
the time-realignment process as described in Section HI of the report. Our 
main program to do this is (arbitrarily) divided into two major subroutines, 
PARTI and PART2. PARTI determines the instantaneous tracker error estimates 
at the 104 times per frame for the S&H data, and PART2 is Intended to make 
use of those time-error estimates. 

For now, PART2 only prints out the 104 per frame time corrections. The 
phrase ^vertical addition” has been used in the past to denote the summing 
into time bins of the properly time-realigned sets of S&H data, and a sub- 
routine VTADD (with entry points VTADI and VTADO) is provided to accomplish 
this. If the waveform computations were all to be done here, we would modify 
PART2 to incorporate VTADD, but this becomes intimately involved with specific 
details of where the S&H data are stored and how to access those data. In- 
stead, the work has been continued at Wallops Flight Center where a single 
Tpa-ffi program now incorporates the two programs PARTI and PART2 plus the 
--''tun -ar'* .details Xo,.vUse .V-TAD]l, documents PARTI 

and PART2 as of approximately February 1974; since then, the development and 
modification work has been done at Wallops. The present Wallops programs 
differ only slightly however from the material presented in this Appendix. 

The organization of the rest of Appendix D is as f oUows : program de- 

scription, list, and flow chart for PART2; then description, list, and flow 
for POLRG, the polynomial regression routine used by PARTI; source listings 
only for GDATA, MINV, MULTR, and ORDER which are needed by POLRG; a descrip- 
tion, listing, and flow chart for SPLINE, the third-degree spline interpola- 
tion used by PARTI; description, listing, and flow chart for WEIGHT and 
XCNVLV; and finally the general description, listing, and flow chart for 
VTADD, together with a note on timing details. Notice that while PARTI re- 
quires a subroutine FDATA to obtain the altitude data upon which PARTI will 
supply no de taxis of FDATA. This is because FDATA again Is spe 
dependant upon how the altitude data is stored or how it is avail- 
able to the program. It is important to notice. that the altitudes returned 
by FDATA are to be in nanoseconds (i.e., the two-way ranging time). 
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Subroutlne Name and Argument List 
Subroutine PARTI (PR) 


Description 

PARTI is the first of a pair of subroutines (PART2, called from PARTI, 
is the other) performing the individual sample— and— hold (S&H) tracker error 
calculation as of February 1974. The 104 per frame tracker errors, in 
nanoseconds, are printed out; the "vertical averaging” of the individual 
sets of 8 S&H results has not been implemented locally (this has been 
accomplished however at NASA/Wallops, and involves relatively few changes 
and additions within PART2). 


PARTI reads general problem parameters and then acquires altitude 
data, an entire submode at a time, from FDATA. A least-squares, polynomial 
regression analysis determines a set of coefficients COE which characterize 
the fitted altitude over the submode. These coefficients COE and the 8 
per frame altitude averages are used to derive average altitude residuals 

. ^ V , . . c Jf ~ J 4-^ B -Fl-aTTIO 

'l(^S 'per rrsme,*'; -osev-^. •• -r , - . 

instantaneous tracker errors which are used to set up a spline fxmetion. 
PART2 then will use the spline in determining the 104 per frame tracker 
errors. The weighting function is based on the work described in the "Task 
A Final Report", (Contract NAS6-2307, Applied Science Associates, approx. 

July 1973) . 


Subroutine Arguments (and Dimensions) in Calling List 


PR(20) *= One line of problem identification to be printed out 
by a 20A4 Format. 


Input/Output Performed by Subroutine 

A number of problem parameters are read, and a general problem heading 
is written each time through the subroutine. The principal input 
to PARTI comes through the call to FDATA which provides the altitude data 
input for the entire problem. The quantities read in directly in PARTI 
are listed below; 
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ML ■ Lower limit to the degree of fit to be employed. POLRG 
vlll perform e regression analysis for all degrees 
between ML and M, starting at ML and incrementing up to 
M (unless there is no improvement in the sum of errors 
squared for a degree less than M in which case the 
analysis stops at that point) ■ 

M “ Maximum degree for the least— squares polynomial regres 
sion fit to be performed by POLRG, called by PARTI. 
(Reading M=0 will retain the already-used values for 
M,ML,NPLOT,NTABXY, and NTEST.) 

NPLOT =s Integer describing how the results of POLRG are to be 

printed out, according to the rules established by POLRG: 

NPLOT = 5 - Print a table of data input, fitted 

value, and residuals for each separate 
degree of the regression, 

NPLOT = 4 - Same as NPL0T=5 but in addition call a 

. PLOT subroutine at each separate degree, 

NPLOT = 3 - Print the table of values and residuals 
only for the final degree fitted, 

NPLOT = 2 - Same as NPLOT-3 but in addition call 
PLOT at the final degree, 

NPLOT = 1 - No table of residuals, no plot. 

NTABXY » Integer determining whether to print out table of weighted 
vs. unweighted tracker error estimates, 8 per frame, in 
PARTIA: 

NTABXY » 0 - No printout 

NTABXY 0 - Print table 

NTEST » Integer setting how much of the entire problem is done: 

NTEST - 2 - Go from start through POLRG only 

NTEST - 1 - Go through the weighting process (which 
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estimates instantaneous tracker, given 
the 8 per frame average tracker errors) 
NTEST = 0 - Do entire problem including determining 

and printing out the 104 per frame tracker 
errors (in nanoseconds) . 



X(2200) X contains the one-per-frame averages of the 8 per frame 
^ altitudes and also contains the times corresponding to 

these frame-averages. The rules for location of these 
input data to POLRG are set by the requirements of GDATA 
called by POLRG: FDATA puts the data for an entire sub- 

mode into the required form. N is the number of such 
one-frame altitude averages. Once the coefficients COE 
have been determined X becomes ^necessary, and it is used 
as temporary storage for the option NTABXY?^0. 

..T(1600) .EDATA also .. returns, the 8 .per frame altitudes in array Y 

Y(1600) with the corresponding times in array T. NY is the total 

NY number of such points, and will be equal to the number 

of frames times 8. 

COE (11) This array contains the coefficients determined by the 
least— squares polynomial regression; hence, all the 
smoothed altitude history for the entire submode is 
contained in COE upon return from POLRG. 

Other Subroutines Called 

FDATA(M,N,X,T,Y,NY) - Obtains the altitude data, an entire submode 

at a time. M is necessary in the argument 
list in order that FDATA stores the one-per- 
frame altitude averages in X. 

POLRG (N,X, ML, M,NPL0T, COE) - Performs the least-squares polynomial 

regression analysis on the frame-average 
altitudes and times in X, to return the 
polynomial coefficients COE. 
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WEIGHT(Y,NY) — Performs the weighting function necessary 

to convert average tracker errors to estimate 
instantaneous tracker errors. The averages 
in Y are converted to Instantaneous estimates 
In Y upon the return from WEIGHT. (Since 
the Y contents are changed in value by calling 
WEIGHT, the Y values are copied over into X 
before calling WEIGHT, in case a **before-and- 
af:te.r” printout is desired as signalled by 

SPI.INE(T,Y,NY,TI,SI,0) < Calling SPLINE with the last argument zero 

uses T,Y,and NY as the input to set up spline- 
fit coefficients within subroutine SPLINE; 
subsequent calls with the last argument equal 
1 will return in SI the spline result corre— 

/ spending to an input time TI, TI and SI are 
'““tsttused arguments in the spline set-up step. 






• DOHBLE PP7CI3T0V X ,T , Y,C0Z,T:, n ^ 

DPkTA *^L/5/,'^/“/»N?LOT/C/#*!T?.?^XY/C'/,HT*'?T/0/,SGt DO/ 


15 FDRMA?(20T3l 

IF (K.t^.O) GO ^9 


- HTABXY^KTRDXY 
■ MTFST-KTr:ST 


ii To^kt"^ {//* •,2c?u/* psobl::^ pA??inFrz?s: ml-* • ,h=* ,1**# 

1 «*i:pioT=*,i“#’r'i^>PXY=« and k7E?t=* ,ia//) 



ISO FOFtlA?(//» ... !13B1AL ?S7n??l TO njlll PFOG,..V/) 
BETOFH ~ ■ 


' I I — 5TT3^r— T7> , 1 

-■—18 POPKRT(//» .FPPOS'RFTUEN TO RAIN PPtSGi'.V * //) ' 
ffPTTF (3,20) (P” (T) # 1=1 ,20) , 1L» R, SFLCT, aTABXt,VTFST 


200 DO 155 X = i,f1?LS1 
155 COF{T)=O.DO 




-DO 250 1=1, VY 
*SI=COE (RP1S1) 


* A — A \ 9 

DO. 225 J=1,1 

JH=HPLS1-J 


SI=tY{i)-si 

T(I)=ST ' 

2 5 0 X { 11^^ 

IF (NTZST.F3.2) 30 TO 239 
CAtL hztght(Y,j:y> 


— fTmyXT.:4;:« I --xi::. i-,.. 

257 POPf*AT(//* AHZ I:X{I);Y(I) ^ * -*/ 

■1 (• ,^13.6, «;• ,013.6,2X1 n 





- CALL PAPT2 (T,Y,NYI 
239 WPIT** (3,3'»'') 


.V* ‘ » \r ^ 

CO TO 13 

END _ _ 
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Subroutine Name and Argument Lls^ 
Subroutine PART2(T,Y,NY) 


Description 

PART2 continues from PARTI the individual S&H tracker error cal- 
culation (104 per frame) as of Feb, 1974. PART2 uses the spline coeffi- 
cients set up in PARTI; for now, the 104 per frame tracker errors are just 
printed out, but eventually the vertical averaging will have to be per- 
formed over a frame of data. This will necessitate modifying PART2. (Such 
modifications are being carried out at NASA/Wallops.) 

Subroutine Arguments (and Dim e nsions) in the Calling Lis t 

Y(1600) = The 8 per frame instantaneous tracker error estimates, 
are used to set up SPLINE. Y is used in the computa- 
tions within SPLINE. 

T(1600) « The times corresponding to the Y data. 

NY « Number of eieiaents? xii — 

data frames. 

Input /Output Performed by Subroutine 

The 104 S&H tracker error corrections per frame are printed out 
as are the 8 per frame values of Y. the tracker error input data to 
SPLINE. The time at the start of each frame is also printed out. 

Other Subroutines Called 

SPI.INE(T,Y.NY.TI,SI,1) - With the last argument in the list 
non-zero, SPLINE returns in SI the spline-determined tracker error 
corresponding to input time TI. 


Qtrher Important Subroutine Vari ables (and Dimension ^ 

PRNT(8) “ An array used to store 8 successive values of the 
tracker error in order to print one line- There 
will be 13 such lines for each frame, since (13 lines/ 
frame) times (8 values/line) equals 104 values/frame. 
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N ■* Integer from 1 to 104 specifying which of the 104 S&H 
times in a given frame is to be found. Notice that 
the calculation of TI in PART2 makes specific use of 
the truncation which occurs in integer division in 
FORTRAN. 
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-sv r t ; ‘Y / vr i 

DI?*SNSI^^ T (‘iO'') , y (800) f20) ' *■ - • 

000BL3 P^FCrsTOM T,Y,7I,Sr,r?WT,TE,CA,C3 

.-i 


T3=T(1J+0. 976nC 
DO 300 K=2,f:F?l 


DO 3 20 
jsB* (K-. 1 ) fl 


«BITF (3,330) K,73, (PFKT(T) , 

330 FOP«AT(/» FRAi? t* , 13 ,' ,ST'.?T-FEiMF: T= • , F 8 . 3 , • :L*?F 7 CCl=*, 
1 i geti »,-q ? '/ Q (ix^f 7 » 2 h 


DO 350 SD5=1,13 
DO ftOC J=1,9 




HD=N /2 

QA-HD 




TI-TB+.004PO* (0B*3 . D0-pA-2.DC) 

, CALL S?LI«3(r,Y,'iy,TI,SI/1) 

r ^ ^ ■ _ 

• •• - I'M 1 • w ^ 


‘ V y 


H«=H 

SL*HO-^ — - • 

WPITF-O (PP>;T (I) 




360 FO®?.AT(* ' ,T«,8(1X,F7.2)) 

350 CONTI!ri?» 

■J ^O — , 0 6 PQ 


PBTUPN 

FHP - 
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Sub routine Name and Ar^^ument L ist 

POLRG (N , X , ML ,M, NPLOT , COE) 

Description 

POLRG is a double-precision and modified version of the sample program 
of the same name in the IBM Scientific Subroutine Package, Version III, 

(IBM Form H20-0205-3) . This is a least-squares polynomial regression 
subroutine using other IBM Scientific Subroutine Package subroutines to 
perform the analysis. Our POLRG is called from subroutine PARTI of the 
Skylab tracker error program, but may be used for other general polynomial 
regression applications. 

Data is carried into POLRG by the array X, and a regression analysis 
Is performed for progressively increasing degree of polynomial (starting 
from degree ML) up to a maximum degree M. Printout or plotting is control- 
led by NPLOT. If there is no reduction in the residual sum of squares 
between two successive degrees of the polynomials, the subroutine terminates 
the problem before completing the analysis for the highest degree poly- 
nomial specified. When the problem is either completed or terminated, 
the polynomial coefficients are returned in array COE. 

Subroutine Arguments (and Dimensions) in Calling List 

X(1100),N - The input data is supplied in the double precision 
array X. The number of input variable pairs is N. 

The individual values of the Independent variable 
are stored in the first N locations of X, or in 
the first column if one views X as a matrix of dimension 
N by M+1. The Individual values of the dependent 
variable are stored in the last column of this matrix. 

The dimension of X must be equal to or greater than 
N*(M+1) . 

KL - ML Is the degree of polynomial at which the problem 
starts. ML must be equal to or less than M; the 
problem starts at ML and increments up to M(unless 
no improvement In the residual sum of squares is 
reached, terminating the analysis before M is reached). 



-82- 


M The maximum degree specified for the polynomial. 

NPLOT - Integer describing how the results of POLRG are to be 

printed outf according to the rules establxshed by POLRG. 

NPLOT = 5 - Print a table of data input, fitted 

value, and residuals for each separate 
degree of the regression, 

4 - Same as NPL0T=5 but in addition call a 
PLOT subroutine at each separate degree, 

3 - Print the table of values and residuals 
only for the final degree fitted, 

NPLOT = 2 - Same as NPL0T=3 but in addition call 
PLOT at the final degree, 

NPLOT = 1 - No table of residuals, no plot. 


NPLOT = 

REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR 

NPLOT == 


by POLRG. If we let XI be an individual value of the 
Independent variable, and YI be the corresponding value 
on the fitted curve, then 

YI = COE(l) + C0E(2)*XI + C0E(3)*XI**2 

+ * . . + COE(11)*XI**10. 

The dimension of COE must be at least Mt-l. 



The output from POLRG is controlled by NPLOT as described above. 


Other Important Subroutine Variable 

DI(100),D(66), - These are various work spaces. DI,D,B, 

B(10),E(10),SB(10), E,SB,T,XBAR,STD,SUMSQ, and ANS are all 

T(10),XBAR(11) ,STD(11), double-precision variables. The dlrnen- 
SUMSQ(ll) ,1SAVE(11) , slon of DI must be at least M*M. The 

ANS(IO) dimension of D must be at least (M+2)ife(M+l)/2, 

The dimensions of B,E,SB, and T must be at 



-83- 


least M, and the dimensions of SBAR,STD, 
SUMSQ, and ISAVE must be at least M+1. 

Our POLRG is set up for a maximum M-10 
and converts any higher value to 10, and 
for M=10 the maximum N allowable (for 
present dimension) is 100. 

P(600) - P carries information for the plot of 

observed data aiid/or printout of estimates. 
The dimension of P should be at least 
3*N if a plot is degired. 

Other Subroutines Called 

The following subroutines are all from the IBM Scientific 
Subroutine Package, and are unchanged except that they 
have been converted to double precision. The descrip- 
tions are from the Scientific Subroutine Package 
documentation and are repeated here for convenience of 
reference . 

GDATA(N,M,X,XBAR, STD, D, SUMSQ) - This subroutine generates independent 

variables up to the Mth power and calculates means, 
standard deviations, sums of cross-products of deviations 
from means, and product moment correlation coefficients. 

ORDER (MM, D, MM, I, ISAVE, DI,E) - Constructs from a larger matrix of correla- 
tion coefficients a subset matrix of intercorrelations 
among independent variables and a vector of intercor- 
relations of independent variables with dependent 
variable. 

MINV(DI,I,DET,B,I) - Uses standard Gauss- Jordan method to invert a matrix 

and calculate the determinant . 

MULTR(N, I, XBAR, STD, SUMSQ, DI,E, ISAVE, B, SB, T,ANS) - Performs a multiple 

regression analysis for a dependent variable and a set 
of Independent variables. 
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PLOT(LA,P,N, 3,0,1) - A special plot subroutine can be provided for 

the input data and the fitted resulting curve. 
The IBM documentation provides a PLOT subroutine 
but in our modified POLRG it is important to 
notice that P is a double-precision array. 

Most of our waveform work has been done using 
a dummy PLOT routine. 


Remarks 

Following first the source listing and then the flow diagram for POLRG, 
we provide the source listings for GDATA,ORDER,MINV , and MULTRj for any 
further details of these routines, reference to the IBM documentation is 
necessary . 







’~^S0BP0f1?Z9*. 

PI-EVSTOM (11) ,ST5(11) ,CO»(11» ,SUf*S0(11) ,lSXVr(11) 

0T"EVST1» >»I5 (^0) , P (f 00) 

G4_«t-*- 0A P ,S Tft ir»,-f4 »1 g Cr^7 -r > r8r SPfg ,> > i 5 r&Si 

D07317 PHECr^irP** P^PESTO 
5 pr'?r. K? { ' 1 POLY'i'^f’.T R L ^EC.^fSSION OP PPO?*?*/!!) 


7 PO1iS»r(*0 F~:r;P35tOS COSFFICrrNTSV(''020.12n 

• y rop*rTT^*TTJT7n'r^tnrsT yijr vr. PT >i.M. E r c »p' ; 7 u - ; — — 

1 • rsT°PE r''iY*i7.'‘r A LV) 

9 POFP >•.7 |» 0 • , 'X , * >'?31CP OP ViPD^TIOH* ,7X,' DFGPPE OP»,7X,*S«H 0** , " 
p^-v^r/^iTpr^'jrrrT^rrfTPTMS'/S'x; »rpt* top ' , 

2 B*,*!>j-JX?.rF*,'?x,'5yt'^.r’'*,7x, 'VAL 0 P',e»,*op sor of squxpes») 

19 9'«7 T7 :i 3T F'!55IO;; • , ’ 2 X , 76, F17 ,S, F1 U ,S , FI 1 , 5,P20 . 5) 

n ~Toyry^r — rrrtnrrjT— r3?ur' ^ SGpr.ssTOF • 5 ,ph» . 9) ^ 

12 POPIXT (Bl,*TnTXL* ,19X,I6,r17.5///) 

13 F3F"Xr('0 SO T^PPOVS-ESr*) 

iiT-T3*'7X77*~»~v/27r7*'r!r9t?~OF~rr5-im)it LS V/ * OBSt P VX Tictt -wo . * ; «iy; 

1 *F VM-0P» ,7X, 'Y VALU*',7X,'Y ESUflATE ' , 7X, • PESIDUAt V) 
15fO‘-A?(» ',3X,T6,F13,5,FH*.5,F17.5,F1S.5) 


F (7.GT,1C) N = n 

ML C.DATA (N, 3, X,X3AP, 370,0,50.130) 


IP (.rL.LT,2) 00 TO 1C5 
90 110 T«1 ,PL 


(r.L.OT, P) El»3 


t$xvt(i)=r 

CALL O^nrs (11,0,51, 1, ISIV!, 01, E) • 

OTL— 5-"5-r(?:,l,!)5r,3,7J ■ • 

• TAIL -,'tLTF (‘r,I,XBAF,S?P,S01SC*DI,E,ISAVE,8,SB,T,APS) 
UPTTF 13,*^) I 

TT-'TTnnn — i ac;T7 7 ; ' T 70 "" 

130 SII«T?=A59 (4)-*S01 ' " ■ 

!? <5l'i;p) 1<i9,1tt0,150 

Ttti»-y?T7 E ( ^, 1 . 1 ) 

xrXTT=t 

159 VSTT2 (1,6) A'JS (1) 

If riT •— T *“,7T — rT( Jr; - "3 - ^tyT ) 

VPTTE (3,fl) I 

V«TTE (1,9) . 


BFTTF (3,V'> t,1NS(a) ,A»S (6) ,AKS(10),SnBIP 
STsA?*? (1) 


UFITF (3,12) NT,SI?1S0{.":.1) 
CO’(1) sANS (1) 


160 CO*'(Jf1) »B(J} 

LA«T 

510 IP (I--.) 520,5?o'630 


220 KP3*‘!»S 

>'0 230 K-^,B 

- ■■■• ■■ s P 

P (HP2) »COE (1) 


P(5P3) ■P(KP3)+E (t) •COE(JI-I) - 
L«l+P • — - — 


L«W*1 

DO 2U0 K*«1,W 


-Pi H? ) - ) {■(!; ) 

WPTTP (3,1«) 

-HP-^W+1 ' — 

DO 250 K»1,H . • -- 

NP2 = }iP2f1 — 



HFS1D=P (NP2)-P (NP3) 

HPITE (3,15) K,P(K) ,P(NP2),P(MP3),PESID 
p ..„.(,_|,H> 10T T *e >i 2 ) ■ O R- r EQ i a ) ) CM 

IP (KSXIT) 100,209, 100 

CONTINUE 
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GT>ATA (N , X, YBAP'^S'? 

p A P ( 1 - ) , T n ( 1 > - r 

- - DOtIBLE P 8 EC IS TOM X , X B AP ^ STD , D , 

I? (M-1) 105,105,90 , 

9 0 ^ 

DO 100 1=2, M 
l1=L^fM 

DO 1C^ J«1yW 

I. = L1+J - 

too X(I.)~X(K) * X(v7 f- 

105 MM=ri+1 
DF = N 

^ ^ — 

- DO -115 T = 1,MM - 

XBAP. (I) =0, DO 

DO 110 J-1 , N 

L=I+1 

110 XBAP (T) =XBA? (I) +X(L) 

t 4S-XBAF - (T> =XPf . ? 

DO 130 1 = 1 ,n?i - — 

130 STD (I) =0. DP — 

h=-i^*t n \ 1 ) ^nr.) /? — - 

DO 1 50 1=1 , L 
150 D(I) =0, DO 

DO 170 K-1,M ^ r 

L=0 

DU I 7 0 

■ t2=N* (J-1)'+K ^ I 

T2=X (L2) -XBAR (J) 

S TD ( J) = ST D C » r n- 'T2^ 

DO 170 1=1, J 
L1=W*(I-1)+K 

r l=X ( L I ) -XCA P CH 

■' 

170 D (L) =D (L) +T1*T2 
— 

DO 175 J=1,KM 
DO 175 1=1, J 

L=t4- 1 

175 D (L) =D (L) -STO (I) *3TD ( J) /DF 

— cro — ti }Q T =i , Kn 

L=L+I 

STJMSOd) =D (L) 

I flO -- b r D(l) =D5 ^ PT (UAB5 (P 

L=0 

DO 190 J = 1 ,M?1 

■ DO 1 90 T = 1 ,J ~ 

L=L+1 

190 0 <L) =D (L) /(STD (I) ♦STD (.1) ) 

DT-SQT x T (Df- 1 . ) 

DO 200 I = 1,M'1 _ ’ 

200 STD (I) =STD(T) /OF 

—rPTtTp^ : — “ 

END 


KEPRODUCIBTEITY'Dr THE 
ORIGINAL PAGE IS POOR 
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Subroutine Name and Argument List 
Subroutine SPLINE (X,Y,N,TJ,SSJ,MODE) 


General Background 

Given a set of n points (x,y,) , (x^,y 2 ) arranged in ascending 
order in the Independent variable x, a third degree spline is a means of 
Interpolation to find the value of y corresponding to a specified for 
X, < X < X . The program accomplishes this by connecting each pair of 
adjacent points with a section of a third degree polynomial. There will be 
a different polynomial for each interval and these polynomials are "matched 
up" at the (x,y) points by requiring that the first and second derivatives be 
continuous at those points. The spline function is in some senses the "smooth- 
est" Interpolating function. As in all interpolating procedures, the interpo- 
lating function goes through each input point, but the spline exhibits con- 
siderably less oscillatory behavior between input points than does Lagrangian 
interpolation, for example. 


The spline function of this subroutine is derived directly from the 
chapter, "Spline Functions, Interpolacion, and Numerical suaarature, oy 
T.N.E. Greville, in Mathematical Methods for Digital Computers , Vol « ll.» 

A. Ralston and H. S. Wilf, editors, (New York; J. Wiley & Sons, 1967), 
pp 156-168. Our spline differs from Greville’ s in trading running time for 
storage space; since we were concerned about the overall size of our time- 
realignment programs, we recompute each time into the routine several quanti- 
ties which Greville *s program stores in arrays. The index search procedure 
is modified also; since our use of SPLINE will usually be for a series of 
monotonically increasing independent variables, we let the index search start 
at the last index found by the program. 


SPLINE operates in two modes: an Initial settingup with the input points 

x,y; and the return of a value SSJ for an input value TJ of the independent 
variable (corresponding to X). We use an indicator, MODE, in the argument 
list of SPLINE to distinguish between these modes. 


Subroutine Argument (and Dimension) in Calling Llsjt 

X(800),Y(800) = The N input points (X,Y) for setting up SPLINE 
N When called with MODE « 0. The successive values 
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of X must be distinct and must be In ascending 
order. 

TJ,SSJ » Output of value SSJ is returned for input independent 

variable value TJ, when MODE =1. 

Input/Output Performed by Subroutine 

None if TJ (in M0DE=1) satisfies X(l) £ TJ £ X(N) . 

If not, an error message is written indicating that 
TJ was outside the allowed range. 
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— ^snP«6(JTTNr~5*PLTVr: (X ,T HnpF) 

i»Tr7r5-To:r'X'rt^fr7t*t^^<^r7^^r»^‘nT^^^r « « 

._. — P1»FCISTr"I X , T J* XI * XTPI , Y i Y I SSJ, S? »S21, 5- - P *BI, » 

1 oi*5G%. p^LTi »t>ft.s cs 

^-m-TPS A-yri‘ l : r./i 
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IF (y.L?.P>"?) r=0 T3 2 • ■■- 
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3 IIIN1=X(2)-X{1) 

. DE1YIM= (Y(2)-T(1) I /HItll ■ — 


HI-X (T+1)-X(I) 

pri.YI= (Y CT41) -Y{t) ) /HI 

‘-y2 i=2>7 T*^p7rr TImr.r7^ftN-^^*7m 

S2 (I)=S2I 

C(T) =1.‘^r)0*S2T 
— DELI T.^-WELTI 

52 ni?i=ni J' * 

S2(1)=''.D^ 

S2 ( N) - r> « DO — — 

HI?R=1 

5 ETH=0. 

^ — » Uf>i-x t 2)- x { ir — 

6 DO 10 1=2, Nl J 

' HI=X (T+1)-X (I) 

S 1-, 5?2 *- Hl * 1 / (LT *1 f ! tl? " 

S2I=S2(I) 

-7 0 ?* y 2- (I -*) f (SI - .I‘ PO) *S2~(T^ 

8 ABW=:DX3S (W) ■" 

• IF GT.'Ti) ErX = ABW 

1 0 32 (I) -021^^ “ ~ 

13 IP (E"*A. tE, EPSLH) GO TO 10 
IF f”!?® --T^f TO ®95 


## 
W ' 


' no TO 5 

— p — oo h > — I 

--soil FOPKAtI/* HU-.3F?, 3F IT?PATIOyS=* ,I3 ,* ->B D SPSI IOH^ * , H 12 > 5//) 

— m i=Y ■■ - - • 

ItAS ^9 ^ ^ 

PFTOPW 

750 IF (TJ-X(I)) 110,120,^30 

1 1DHT-P — (TJ X (1 H — 1^* **■ , 1 0 , 1 f P ^ ■ 

'"iJSs "oP»lTt*’lFr,!l-isT TJ=*,012.5,»,- ISSS THAS X(D') 

^ OO 

120 SSJ=Y(T1 

P2TU®N 
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— 200 I-Ifl ' • • 1 J 

IP (TJ-X(T)) 2 10,120 ,200 

170 IF (I.’Q.ILAST) 7D TO 173 

XI=X(I) 

{ T H-1 ^ H 

S2I-S2(T) ■ ; 

S2TPl = *^2 (T+1) ■ ■ 

' ^ — '■ ~~ 
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g ubroutlne Namd and Argument List 
Subroutine WEIGHT (Y,NY) 

^neral Description 

WEIGHT is called to apply a weighting function to a set of NY input 
data point. Y. Upon return from WEIGHT, the original Y will be replaced 
by the set of values from applying the NW individual weights to the input Y. 

The subroutine XCNVLV vrLth an entry point XCFILL is used for the convolu- 
tional weighting. Since generally NW will be considerably smaller than NY, 
the total storage needed is only NW+NY rather thaft the 2*NY which would be 
needed for a straightforward digital convolution. 

The number of weights, NW, and the weight values, W, are set by the 
DATA statement in WEIGHT. Note that with NW » 17, if W(9) is set to 1.0 
and the other weights to 0. , the original Y is returned from W. 

Subroutine Arguments (and Dimensions) in Calling List 

Y(800) = Input points (number - NY) entering WEIGHT; on 
return, Y is replaced by the weighted result. 

Other Important Subroutine Variables 

WEIGHTjNW = WEIGHT (NW) contains the weight as set by DATA 
statement. 

TEMP ® A working temporary array. The dimensions of W 

and TEMP need to be at least as large as NW. 

Input /Output Performed by Subroutine 

None 

Other Subroutines Called 

XCNVLV (with entry XCFILL) - Carries out digital convolution when 
called repeatedly with successive points. XCNVLV 
listing only is provided here. 



SOPPOUTIN?) WEIGHT (Y, NY) 


DOUBLE PPECISION y, yi,w,th«p,yout . . . 

DATA NW/17/,W/-.OOlDO,-.en'^SDO,-.OUiSDr,‘-.03DO,-.OUBSDO,-.057DO 


rAKbViTftl 




• •03n0,“,0^45D0,-. 005BDO,- .001D0/ 
IF (KY.GE.NW) GO TO 5 


3 FOnWAT(//’ NO WEIGHTING DONE, BECAUSE NY=SI3,*rBU? 
BETUPN 


DO 10 I=T,NW\>1 
YI = Y(T) 


DO 20 X=NW,NY 
YT^Y (T) 


JJ=I-NWH 
20 Y(JJ) =YOUT 


DIMENSION TEMP (NH) (NH) 

DOUBLE PBECISTON TEMP, H, YT,XT 


1 J. — W m V J ' 

JP=NII+1 

NP=MOD (JXY-1 ,NH) 


TEMP (NH) =XJ 
DO 100 1^1, NH 


100 YI=YH-?EM? (T) ♦H(M) 
PFTU^N 


K=NH-NP 
DO 210 1=1 ,K 


210 YI = YT + '^E«P (.7) *n (M) 


DO 220 1=1 ,NP 


220 YI = YI-HTEMP (T) (M) 

PSTtlPN 


NP=MOD (JXY~1 ,NN) 

IF (NP.GT.^) GO TO 300 


RETURN 

300 TEMPtN^)=^" 


S 

















Yi-vd) 


CUX XCS\XV(YOUT,Yl.I» 
TDff.U.KM) 

Y(JJ)*YOUT 





1 1 j?4.jn‘+i 1 





t '' — 

1 J^I+K 

I bUJ?-I 

1 YUYI+TE>a»(l)*H(K) 

i 

'— -C 

c 



i> 


RETUWJ 
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Sitbroutlne Name and Argument List 

Subroutine VTADI (MB , EPS , SANDM, SUMl , SUM2 , NSUM) 

Entry VTADD(EPS,SANDH) 

Entry VTADO (SUMl, SUM2, NSUM) 

Description 

This subroutine accomplishes the "vertical averaging” of a number 
of sets of sample-and-hold results once each individual sample-and-hold 
time tracker error correction EPS has been found. ^ VTADI initializes 
the storage locations and sets up several interval variables, VTADD 
is used to enter each successive sample-and-hold set SANDH in bins 
appropriate to the individual tracker error EPS, and VTADO gets ready 
for outputting the means, standard deviations, and total number of 
entries at each storage bin. 

Subroutine Arguments (and Dimensions) in Calling List 

.. _ — -» ..4J fcV <->'.• r\P 

The individual bin width is GS/MB (where GS is the S&H 
gate separation, 10 nanoseconds in the case of SKYI*AB 
S-193) . KMB^IO 

EPS — Individual tracker error in nanoseconds at time 
corresponding to the input SANDH. 

SANDH(8) = Input set of 8 sample-and-hold gate readings. 

SUMl (151) = The set of double-precision bins into which the differ- 
ent Inputs are summed; the call to VTADI zeroes all 
these and the call to VTADO converts all the subtotals 
to Individual bin means. 

SUM2(151) = A set of double-precision bins into which the squares 
of SANDH are added. VTADI zeroes all SUM2, and the 
call to VTADO converts the sum of squares into individual 
bln standard deviation estimates, 

NSUM(151) - An integer array accompanying SUMl and SUM2 which 

records the number of entries to each individual bin. 

The call to VTADI zeroes all NSUM. 
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Other Important Subroutine Variables 

GS “ Sample-and-hold gate separation in nanoseconoo 

S-193 Mode 5, submodes 1 and 2, and GS=25.0 for the other 
submodes of Modes 1 and 5) . 

NBB = An integer locating the base bin, calculated from WB in 

VTADI. All 8 S&H locations are calculated relative to_.NBB 
via JB in VTADD. 

Input/Output Performed by Subroutine 

None except through calling list. ‘ 

Other Subroutines Called 


None 


la 


20 


40 

35 

30 


SUBROUTINE VTADI (MB> EPS, SANDH, SUN 1 SUM) 

DIMEMSI CM SAMDH C 8 ) , SUM 1(151), SUM2 C 1 5 1 ) ,N SUM (151) 
DOUBLE PEECiSIOr-J SUM 1 , STJM2, SHI , DSG, XMM 


DATA GS/10./ 

IF CMB.GT. 10) MB=1 0 

NBB=4^fMB+l 

VB=GS/MB 

DO 10 1=1,151 

SUMKI ) = 0,D0 

SUM2(I )=0.D0 

NSUM(I)=0 

RETURN 

ENTRY VTADD(EPS, SANDH) 
Q=EPS/UB-.5 
IF (EPS.GE.O.) Q=G+K 
JB=NBB+IMT( Q) 

IF (JB.LT. 1 ) JE=1 
IF (JB.GT.81) JB=81 
DO 20 1=1,8 
SHI=SANDH(I) 

SUMl (JB)=SUM1 (JB)+SHI 
SUM2(JB)=SUM2( JB)+SHI*SHI 



NSUM( JB)=MSUM( J3)+ 1 


<JB=JB+MB 

RETURN 

ENTRY VTADO(SUMl, SUM2,NSUM) 

DO 30 1=1, 151 

XNS=NSUM(I) 

IF (XNS.LT,!.) GO TO 30 
XMN= SUM 1(1) /XM S 


SUMl(I)=XMN 

IF (XNS.GT.U) GO TO 40 
SG=16.^^SUM2(I ) 

GO TO 35 

DSG=SUM2( I )-XNS^XMN*XMN 
SG=DSG/(XNS- 1 . ) 

SUM2< I )=SQRT(A3S( SG) ) 


CONTINUE • 

retur:^! 

END 
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ester 

WTAM. 


1-EPS/VR-.5 


MB.CT.IO 
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Tlmc and Position Details for Bins in VTADD 

Tlic accompanying sketch shows the bin numbering and positioning 
details. Notice that this is specific to eight saraple-and-hold gates 
08 In the Skylab S-193 experiment (although the changes necessary to 
go to a different number of S&H gates would be. simple ones), and that 
on additional four gate positions are. located on e^ither side of the 
original eight to allow for tracker jitter of up to + 4 GS. 

The total number of ’’vertical averaging" bins is (15*MB)+1, and 
VTADD allows for MB ^10, so that the maximum number of bins, at MB=10, 

Is 151. 

If one uses JSO as an index to denote the original S&H gate numbers 
(JSO-1, . . . ,8) , and JMB as the index to denote the bin number corresponding 
to JSO for a specified MB, the relation between JSO and JMB is 
JMB==(JS0+3)*MB+1 . 

Also, relative to a local time origin t=0 located at the position of the 
first of the original gate positions, the time at the center of any 
specified bin JMB is given by 

tlme=(JMB-4*MB-l)*GS/MB . 



Added 4 poslcloas ' 
to allow for 
tracker jitter 



Added 4 positions 
to allow for 
tracker jitter 


1 




II 


a. 


WB » GS/1. 


10, 




2V. 


z8. 


ft ^ 


25| I 25j 32 


WB = GS/4., 


t «0 
location 


Sketch of Time and Position Relations ariong VTADD Bins 
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